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^ ' Abstract 

. We analyze the representation of A n as a linear combination of 

', A J , < j < k — 1, where A is a k x k matrix. We obtain a first 

*n ■ order asymptotic approximation of A n as n — > oo, without imposing 

any special conditions on A. 
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1 Introduction 

In a recent article pQ, Abu-Saris and Ahmad showed how to compute the 
powers of a matrix without having to compute its eigenvalues. Their main 
result was: 
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Theorem 1 If A is a k x k matrix with characteristic polynomial 

k-l 

P{x) = x k + J2 a i xj > (!) 

3=0 

then, 

k-l 

A n = Y^bj(n)A j , n>k (2) 

j=0 

where 

bj(k) = -aj, < j < k - 1, 6_i (n) = 0, n > k, (3) 
bj(n + 1) = bj-i(n) — ajbk-i{n), n > k, < j < k — 1. 

The purpose of this paper is to find an asymptotic representation for the 
numbers bj(n) as n — > oo, which using (J2J) will give an asymptotic represen- 
tation of A n for large n. Since the coefficients bj(n) depend only on P(x), 
our estimates will be valid for similar matrices. 

The asymptotic behavior of powers of matrices has been considered before 
by other authors. In [3] and Gautschi computed upper bounds for A n 
and ||^4 n || , where ||A|| is a norm of A. Estimates of ||^4 n || were also studied 
in 0, 0, HQ, M and 

In Friedland and Schneider considered the matrix 

B (m) = A m (l+... + A"- 1 ) , m > 1 

where A is a nonnegative matrix and q is a certain positive integer. They 
proved a theorem on the growth of under the assumption that the 

spectral radius of A is equal to one. Powers of nonnegative matrices were also 
analyzed by Lindqvist in jSJ. Rothblum ^U], obtained Cesaro asymptotic 

N 

expansions of where A is a complex matrix with spectral radius less 

than or equal to one. 

This paper is organized as follows: In Section |21 we find an integral rep- 
resentation for the exponential generating function Gj(z) of the coefficients 
bj(n). We obtain exact formulas for Gj(z) and bj(n) in the special case of the 
matrix A having k distinct eigenvalues. We conclude the section with some 
examples. 
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In Section 01 we give an exact representation and a first order asymptotic 
approximation for bj(n), as n — > oo. We consider the cases of simple and 
multiple eigenvalues. Our formulas are relatively easy to implement and offer 
very accurate estimates of bj(n), and therefore of A n , for large n. We present 
some examples for different cases of P(x). 

2 Generating function 

In this section we shall find and exponential generating function for the 
coefficients bj(n). First, let us define the spectral radius p(A) of the matrix 
A by 



p(A) =max{|A| | P(A) = 0} . 



(4) 



Theorem 2 Let Gj(z) be defined by 




(5) 



n>0 



Then, we have 



c+ioo 





c— 100 



where c > p(A) 




0<j<k-l, 



(7) 



1=0 



and P(s) is the characteristic polynomial of A defined in |7p. 



Proof. If we use in (J3J), we obtain 



G' 




Gj(0) = -a j} 0<j<k-l 



(8) 



with = 0. Taking the Laplace transform of Gj(z), 



oo 
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in (jSJ) we get 

sLj + a,j = Lj-i — (ijLk-i, < j < k — 1 (9) 

and L_i(z) = 0. ■ 

The solution of (0) is given by 



- 3 -lPj( S 

^3\ J ) — ° 

where 



Lj (s) = -s k ~^f^, 0<j<k-l (10) 



Pj{s) = ^a l x\ 0<j<k-l, 

1=0 

and P-i(s) = 0. Inverting the Laplace transform in (fTU|) the theorem follows. 
Remark 3 Since 

limL i (s) = 0, < j < A; — 1 

we can replace the Bromwich contour in (OJ) with a circle C of radius R 
centered at the origin 

c 

with R > p(A). 
Corollary 4 If 

P(s) = ( a -Ai)(s-A 2 )...(a-A fc ), 
where the eigenvalues A, are all distinct, then 

Gj(z) = -£ W-'- 1 |?XT ex P ^ ( 12 ) 

1=1 \ l > 



and 

k 



i=i 
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Proof. Applying the residue theorem to ([lip we obtain 



G s (z) 



- £ Res 

P(A)=0 



-lPj(' 



P(s] 



-e zs :X 



which in turn gives (|12|) after computing 



Re s 



P(s] 



-e zs :\. 



s^\t F3K J P(s) 



Writing (0 as 

= -£ (AO 



and changing the order of summation, we have 



n>0 



z 



which implies 



6, (n + *) = -£ M 



z=i 



n+k—j—l 



(14) 



Example 5 In ^ the authors considered the following examples: 
1. 

p(x) = x 3 - 7x 2 + lQx - 12 = (x - 2) 2 (x - 3) . 
[/smg < f7^) we /jai>e 

-s 2 (12) 



^o(s) 

Z^s) 



(s - 2) 2 (5 - 3) 
-s (16s- 12) 
(5 - 2) 2 (5 - 3) 
- (-7x + 16x - 12) 
(s-2) 2 ( S -3) 
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and inverting we obtain 



G (z) = -12 (8 + Az) e 2z + 108e 3z 
G x {z) = 4 (23 + lOz) e 2z - 108e 3z 
G 2 {z) = -4(5 + 2z) e 2z + 27e 3z . 



Expanding in series we get 



G (z) = -96V T— - 48V2 n - 1 rA + 108^3^ 

n>0 n>0 n>0 

and from |pj) we conclude that 

b (n + 3) = -96 x 2" - 48 x 2 n ~V + 108 x 3" 

or 

b ( n ) = -96 x 2 n ~ 3 - 48 x 2 n ~ 4 (n - 3) + 108 x 3 n ~ 3 
= -3(l + n) x 2 n + 4 x 3 n . 

Similar calculations give 

6i (n) =^njx2 n -4x 3" 
(n) = - f 1 + in) x 2" + 3", 
in agreement with the results shown in £7]/. 

2. 

P(x) = x 3 — 5x 2 + 6x = x (x — 2) (x — 3) . 
We can apply 173)) directly and obtain 

b (n) = 

6j (n) = - x 2" - - x 3" 
v ; 2 3 

b 2 (n) = -- x 2 n + - x 3". 
v ; 2 3 
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3. 

P{x) = x 5 -5x 4 +10x 3 -20x 2 -15x-4 = (x - 4) (x 4 - x 3 + 6x 2 + 4x + l) . 

Although (as the authors noted) MAPLE is unable to compute the zeros 
of P(x) exactly, it can provide us with very accurate numerical approx- 
imations 

Ai = 4 

A 2 = 0.8090169944 + 2.489898285i 
A 3 = 0.8090169944 - 2.489898285i 
A 4 = -0.3090169944 + 0.2245139883i 
A 5 = -0.3090169944 - 0.2245139883i 

which we can use in \13]) to get 



b (n) = _ y (A ) n- i-i PjW) 

3K 1 305 l) P'{\ x ) 



1=2 

with 

C = 1,(7! = 4,C 2 = 6,C 3 = -1,C 4 = 1. 
Note that, for < j < 4, we have 

E #H = O (|A 2 D = O (2.618- 



Z=2 

as n — > oo. 



Example 6 Let A be the matrix 

A = 

with characteristic polynomial 



1 2 
-1 -1 



P(x) = x 2 + 1 = (x - i)(x + i). 

t/szng 17^) we /iai>e 

6 (n) = cos (^nj 
bi(n) = sin (77^ 
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and from (dj) we get 
A n -- 

In particular, we have 



cos (|n) + sin (| nj 2 sin 

— sin (|n) cos — sin (|n) 



I, 


n = 0(4) 


A, 


n = 1(4) 




n = 2(4) 


-A, 


n = 3(4) 



where I denotes the identity matrix. 

Example 7 This example appeared in [9]. Let A be the matrix 

Hi 

with characteristic polynomial 

P(x) = x 2 — x — 1 = (x — a)(x — P), 

where 

a = i(i + VB), /J = i(i-VS) 

TTien, /rom 17^) . we /lave 

6 (n) = ^=(a"- 1 -^- 1 ) =/„_! 

bl (n) = ^=(a n -(3 n ) = fn 
where f n is the nth Fibonacci number. Thus, 
A n -- 



fn + fn-1 fn \ __ ( fn+1 fn \ 
fn fn—l) V fn fn— 1/ 
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3 Asymptotic analysis 

We begin by finding an integral representation of the coefficients bj(n). 
Lemma 8 The numbers bj(n) can be represented as 



P{s) 



where C is a circle of radius R > p{A) centered at the origin and the poly- 
nomials Pj{s) were defined in |7p. 



Proof. Since the power series 



n>0 



converges uniformly on \s\ < R, we can interchange integration and summa- 
tion in (fTTj) and obtain 



n>0 



1 

2tt[ 



P s) 



Z" 

n\ 



Then, (0 implies 



n+fc-j-lPj( S ) 



bj (k + n ) = --L[ s ^-^^ds. (J < j < /,■ I 
2-kiJ P{s) 



and the result follows. ■ 

Remark 9 An alternative method for approximating the coefficients bj (n) is 
to write J73J) as 



bj(n) 



R n-j 

2vr 



2tt 



exp [it(n — j) 



Pj(Rt 



P{Re l ) 



-rdt 



(16) 



with R > p(A) and to compute the integral Mb]) numerically. This approach 
offers the advantage of avoiding the computation of the eigenvalues of A. 
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We now have all the necessary elements to establish our main theorem. 
Theorem 10 Let 

p(A) = |A| > |A 2 | > • • ■ > |Ar| 
be the eigenvalues of the matrix A, i.e., 

P{x) = (x - A) m (x - A 2 ) ma • • - - X r ) mr (17) 

with r < k. Then, 

6 ( n ) ^ _ A n-m-, JiW Jn-k\ n ^ 

n ; P( m )(A) \m-lj K ' 

•where 

Am p 

PW(A) 



Proof. To find an asymptotic approximation of (fT3j) . we shall use a modified 
version of Darboux's Method 1121. We write 



P(s) P(s) 



so that deg (s k J 1 pj) = k — 1 and deg (P) = A;. 
From (fTTj) we have 



P( S ) (s-X) m g(Xy 
where 

^(x) = (a; - A 2 ) m2 • • • {x - A r ) 
Using the Binomial Theorem, we obtain 

71— fe 

( .s- — \ )' ( " / 



a"-* = £ (a - A) 1 ( " ^A""*"'. (20) 



Combining (fTTJJ) and (f2TIJ) . we get 

„n-j-lPi( g ) _ A n-m-j M^l f w ~ * \ s , A 



P(s) ~ ^(A) U - V - A) 
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and therefore 

J 27ri7 P(s) g(X) \m-lj' K ' 

c 

To find the value of g(X), we use L'Hopital's Theorem 

P(s) p( m Us) P^(X) 

g(X) = lim^L = limi— ^ = f— fl. (22) 
s^\[s — X) m\ ml 

Replacing (J22I) hi (jZED, we obtain (HHJ) . ■ 

Remark 11 Note that when m = 1 we recover the leading term in M'J\) . 

If more than one eigenvalue has absolute value equal to the spectral ra- 
dius of A, the asymptotic behavior of bj(n) can be obtained by adding the 
contributions from each eigenvalue. We state this formally in the following 
corollary. 

Corollary 12 If 

p(A) = |Ai| = |A 2 | = ■ • ■ = |A r | , 
with respective multiplicities mi,m2, . . . , m r , then 

^^-pX^^^i^.^), (23) 

Remark 13 Since 



n — k 



mi - lj (mi - 1)! 



n mi \ n — > oo 



we have 

b 3 {n) ~ -J>)— j^^rn^-\ n-oo 



i=i 



Therefore, in the case of several eigenvalues located on the circle \s\ = p(A), 
the dominant term in \2ty) will correspond to the eigenvalue with the greatest 
multiplicity. 
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Example 14 In Example^ (1) we consider 

p{x) = x 3 - 7x 2 + 16x - 12 = (x - 2) 2 (x - 3) . 

In t/iis case, A = 2, m = 2 and A; = 3. From ilty) . we get 

bo(n) ~4x3 n 
b^n) ~-4x3 n 
62 (n) ~ 3 n 

which are the leading terms in the solution previously obtained. 

Example 15 We now consider the case of more than one eigenvalue having 
absolute value equal to p(A). Let 

P(x) = X A + x 3 - 15x 2 - 9x + 54 = (x - 2) (x - 3) (x + 3) 2 . 

In this case, Ai = —3, m\ = 2, A2 = 3, m 2 = 1 and k — 4. From we 
nave 

6d(n)~~n(-3) B + |(-3r-i3- 

61 (n) ~ — ra (-3) n - - (-3) n - — 3 n 
u ; 10 V ; 5 V ; 12 

b 2 (n) ~ — n (-3) n - — (-3) n + -T 
v 7 45 V ; 45 V ; 9 

63(71) ~ -— n (-3)" + — (-3) n + —3". 

3V J go V ; -r 45 V ; -r 3 g 

The exact values are 

b (n) = --n (-3) n + — (-3) n - -3 n + — 2 n 
uv ; 5 V ; 50 V ; 2 25 

&!(„) = — n (-3) n - — (-3)" - —T + — 2 n 
u ; 10 V ; 300 V ; 12 25 

62 (n) = — n (-3) n + — (-3) n + -3 n - — 2 n 
v ; 45 V ' 225 V 7 9 25 

63 (n) = -— n (-3)" + — (-3)™ + — 3 n - — 2 n . 
3V ; 90 V ; 900 ^ ' 36 25 

As we observed before, the main contribution comes from the eigenvalue of 
maximum multiplicity, in this case Ai = —3. 
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Example 16 Finally, let's consider the case of complex eigenvalues of mul- 
tiplicity greater than one located on the circle \s\ = p{A). Let 

p{x) = x 5 - 9x 4 + 34x 3 - mx 2 + 65x - 25 
= (x-l)[x-(2 + i)] 2 [x-(2-i)] 2 . 

In this case, \\ = 2 + i, mi = 2, A2 = 2 — i, m>2 = 2 and k = 5. From \2ty) . 
we obtain 



b (n 
bi(n 
b 2 (n 
b 3 (n 
biin 



- (v^j {n — 5) [cos (On) — 7 sin (On)] 
-— (y^Y (n - 5) [2 cos (On) - 39 sin (On)] 
-— (VEy '(n-5) [2 cos (0n) + 31 sin (0n)] 
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a/5) (n-5)[2cos(6>n) + llsin(#n)] 



~ (V5^ n (n-5) [cos (0n) + 3 sin (On)] 



with 



= arctan - 1 . 

2 



The exact values are 

b (n) = ^ (ytj n [(n - 21) cos (On) + (-In + 22) sin (On)] + 



b x (n 
b 2 (n 
b 3 (n 
64 (n 



25 
T 



-— (ytj n [2(n - 50) cos (On) + (-39n + 125) sin (On)] - 10 
-— (v 7 ^)" [(2n + 65) cos (On) + (31n - 100) sin (On)] + y 
— (ytj n [2(n + 10) cos (On) + (lln - 35) sin (0n)] - 2 
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